%**************************************************************************
%
%           "Endogenous Liquidity and Capital Reallocation"
%      Cui, W., Wright, R., & Zhu, Y. (2024), Journal of Political Economy
%                       Last Modified: Feb 2024.  
%
%**************************************************************************

%%%% Run business-cycle simulations for the main (CWZ) model

clc; clear; close all; format short



%% Set parameters related to the steady state
Main_iid_calibration


%% 0. ------------------------control of experirements---------------------
simulate_           = 1;
monetary_shock      = 0; %not in the paper


%% 1. ------------------------- Initilization------------------------------
load param_calibrated
load eqm_calibrated
par = par_calibrated;
eqm = eqm_calibrated;

%%%% store the steady state values
K_bar            = eqm.K_bar_K * eqm.K;
E_a              = eqm.E_a_K * eqm.K;
E_p              = eqm.E_p_K * eqm.K;
R_ratio          = (E_a + E_p) / (E_a + E_p + par.DELTA * eqm.K);
P_share          = E_p / (E_a + E_p);        
par.B_SS         = eqm.B;
par.Z_SS         = eqm.Z;
par.K_SS         = eqm.K;
par.H_SS         = eqm.H;
par.W_SS         = (eqm.B / (1 - par.ETA))^((par.ETA - 1) / par.ETA) * par.ETA;
par.L_SS         = eqm.L;
par.Y_SS         = eqm.Y;
par.C_STAR_SS    = eqm.c_star;
par.K_BAR_SS     = K_bar;
par.ALPHA_0_SS   = eqm.ALPHA_0;
par.q_capital_SS = eqm.q_capital;
par.q_money_SS   = eqm.q_money;
par.q_profit_SS  = eqm.q_profit;
par.surplus_SS   = eqm.surplus_net;
par.q_liq_SS     = eqm.q_liq;
par.E_a_SS       = E_a;
par.E_p_SS       = E_p;
par.R_ratio_SS   = R_ratio;
par.P_share_SS   = P_share;
par.DM_price_SS  = eqm.average_2nd_price;
par.prod_disp_SS = eqm.prod_std;
    
% %%%% Options value to be used in fsolve
% options                     = optimset('Display','off'); 
% options.OptimalityTolerance = 1e-6; 
% options.StepTolerance       = 1e-6;
% options.MaxIter             = 1000000; 
% options.MaxFunEvals         = 1000000; 
% options.FunctionTolerance   = 1e-6;


%% 3. ----------------Run CWZ model -----------------------, dynamics
disp('**********************> Model 6 (CWZ) Simulation <**********************')

%%%% Related to dynamics
par.RHO_CREDIT      = 0.83;
par.RHO_A           = 0.83;
par.RHO_MU          = 0.83;
par.RHO_G           = 0.83;

%%%% run DYNARE file to obtain simulations
if simulate_ == 1

    cd dynamics_iid        
        disp('Simulatino with Credit and Productivity Shocks')
        par.SIGMA_CREDIT    = 8.21/100;
        par.SIGMA_A         = 2.76/100;
        par.SIGMA_MU        = 0;
        par.SIGMA_G         = 0;
        save params par   
        dynare model_6_IRFs.mod nostrict
        oo_credit_prod = oo_;

        disp('Simulatino with only Productivity Shocks')
        par.SIGMA_CREDIT    = 0;
        par.SIGMA_A         = 3.83/100;
        par.SIGMA_MU        = 0;
        par.SIGMA_G         = 0;
        save params par   
        dynare model_6_IRFs.mod nostrict
        oo_prod = oo_;        
    cd ..

    disp('**********************> Results <**********************')

    
    fprintf('output volatility under A&credit is %.4f \n \n', sqrt(oo_credit_prod.var(1)))
    fprintf('output volatility under A is %.4f \n \n', sqrt(oo_prod.var(1)))
    
    std_dev_credit_prod = sqrt (diag(oo_credit_prod.var)) ./ sqrt(oo_credit_prod.var(1)); %the first position is output. Consumption starts from the 2nd
    std_dev_prod = sqrt (diag(oo_prod.var)) ./ sqrt(oo_prod.var(1));

    corr_w_output_credit_prod = oo_credit_prod.contemporaneous_correlation(1,:); 
    corr_w_output_prod = oo_prod.contemporaneous_correlation(1,:); 

    corr_w_pi_credit_prod = oo_credit_prod.contemporaneous_correlation(8, :); % inflation is 8th position
    corr_w_pi_prod = oo_prod.contemporaneous_correlation(8, :); % inflation is 8th position

 
    disp('-----------------------------')
    disp('Table 7 result and more (last three rows, some are discussed in the paper):')
    disp('Relative Std.Dev. and Correlation')
    Table_7 = array2table([std_dev_prod(2:11), std_dev_credit_prod(2:11),...
       corr_w_output_prod(2:11)',  corr_w_output_credit_prod(2:11)'],...
       'RowNames', {'Consumption','Investment','Hours','TFP','R share', 'P share',...
        'Inflation','ALPHA_0','Prod disp','2nd Mkt Price'},...
       'VariableNames', {'std.dev. (A)', 'std.dev. (A&Credit)',...
        'corr (A)', 'corr (A&Credit)'})

    disp('-----------------------------')
    disp('Corr with inflation')
    array2table([corr_w_pi_prod(1:7)', corr_w_pi_credit_prod(1:7)'],...
        'RowNames', {'Output','Consumption','Investment','Hours','TFP',...
        'R share','P share'},...
        'VariableNames',{'corr (A)','corr (A&Credit)'})

    disp('-----------------------------')
    disp('Corr of hours with TFP under two shocks and only productivity shock') % Note the data is 0.48
    [oo_credit_prod.contemporaneous_correlation(5,4), oo_prod.contemporaneous_correlation(5,4)]

    disp('-----------------------------')
    disp('Corr of hours with wage under two shocks and only productivity shock') % Note the number calculated from Christiano and Eichenbaum (1992) is 0.16
    [oo_credit_prod.contemporaneous_correlation(4,12), oo_prod.contemporaneous_correlation(4,12)]
end

if monetary_shock == 1
    %%%% Impulse response for monetary shock
    par.SIGMA_CREDIT    = 0;
    par.SIGMA_A         = 0;
    par.SIGMA_MU        = 0.01; % 1% money growth shock
    par.SIGMA_G         = 0;
    cd dynamics_iid 
        save params par   
        dynare model_6_IRFs.mod nostrict
    cd ..    
    xgrid = [-5:20];
    
    ss_value = zeros(6,1);
    figure()
    subplot(2,3,1)
    plot(xgrid, [ss_value; G_pi_e_Mu] * 100)
    ylabel('%')
    title('Inflation')
    
    subplot(2,3,2)
    plot(xgrid, [ss_value; Mu_e_Mu] * 100)
    ylabel('%')
    title('Money growth')

    subplot(2,3,3)
    plot(xgrid, [ss_value; G_INV_e_Mu] * 100)
    ylabel('% deviation')
    title('Investment')
    
    subplot(2,3,4)
    plot(xgrid, [ss_value; G_output_e_Mu] * 100)
    ylabel('% deviation')
    title('Output')

    subplot(2,3,5)
    plot(xgrid, [ss_value; G_R_ratio_e_Mu] * 100)
    ylabel('%')
    title('R share')

    subplot(2,3,6)
    plot(xgrid, [ss_value; G_P_share_e_Mu] * 100)
    ylabel('%')
    title('P ratio')
end

